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The Sensitivity of the Parkes Pulsar Timing Array to 
Individual Sources of Gravitational Waves 
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ABSTRACT 

We present the sensitivity of the Parkes Pulsar Timing Array to gravitational waves 
emitted by individual super-massive black-hole binary systems in the early phases of 
coalescing at the cores of merged galaxies. Our analysis includes a detailed study of the 
effects of fitting a pulsar timing model to non-white timing residuals. Pulsar timing 
is sensitive at nanoHertz frequencies and hence complementary to LIGO and LISA. 
We place a sky-averaged constraint on the merger rate of nearby (z < 0.6) black-hole 
binaries in the early phases of coalescence with a chirp mass of 10^" Mq of less than 
one merger every seven years. The prospects for future gravitational- wave astronomy 
of this type with the proposed Square Kilometre Array telescope are discussed. 

Key words: gravitational waves - pulsars: general. 



1 INTRODUCTION 

In the era of ground- and space-based gravitational-wave 
(GW) detectors, GW astronomy is becoming increasingly 
important for the wider astronomy and physics communi- 
ties. The ability of the current GW community to provide ei- 
ther limits on, or detections of, GW emission is of enormous 
importance in characterising astrophysical sources of inter- 
est for further investigation. It is possible that GW detection 
will provide the only means to probe some of these sources. 
The sensitivity of existing and future observatories to indi- 
vidual GW sources, such as neutron-star binary systems and 
coalescing black-hole binary systems, has been calculated 
in the ~kHz and ~mHz frequency ranges. The sensitiv- 
ity curves of the Laser Interferom eter Gra vitational- Wave 
Obse rvatory (|Abbott et al.l I2009l f1. Virgo (|Acernese et all 
l2006l rl and the Laser Interferometer Space Antenna (LISA: 
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iLarson et al.ll200Cl lJ cover these frequency ranges. The sen- 
sitivity of GW detectors to individual sources of GWs at 
lower frequencies has not yet been presented in detail. 

Radio observations of pulsars have long b een proposed 
as a means of d e tecting low-frequency GWs llSazhinlll978l; 
'Petweilei' 1 19791: iHellings fc Downs! 1 19831 : IWvithe fc b'oebl 
2003; Jcn et et al.ll2005h . Pulsars are used as a GW detector 
via comparison between a model for their pulse arrival times 
and high precision measurements of these "times-of-arrival" 
(TOAs) at a radio telescope over a period of ye ars (see, e.g., 
iLorimer fc Krameij[20o3 : lEd wards et al.ll2006l 'l. Pulsar tim- 
ing is mo st sensitive to GWs in the ~ nHz fre quency range 
(see, e.g.. |Jenet et al.ll2005l : [Hobbs et al.ll2009l 'l. The sources 
most likely to produce a detectable GW signal in this fre- 
quency range are super-massive black-hole binary systems 
(SMBHBs) in the early phases of coalescence at the cores of 
merged galaxies (see, e.g., Scsana et al. 2 009|). 

Earlier work JRomani fc TavloiJll983l : lKaspi et al.lll994l : 
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iLommenl 12003 : Ijenet et al.1 120061 ) aimed to limit the am- 
plitude of the stochastic background of GWs, either using 
observations of an individual pulsar, or using precise and 
contemporaneous timing of several pulsars (a "pulsar tim- 
ing array" [PTA] ). Here we de scribe methods t o limit or 
detect individual GW sources. In lWen et al.1 (|20ld ). the non- 
detection of GWs from a sing le SMBHB i n the p ulsar tim- 
ing observations presented in iJenet et al.l (|2006l ) was used 
to place limits on the coalescence rate of SMBHBs with 
a range of reds hifts and ch ir p ma sses. The timing residu- 
als reported in Ijenet et al.l (|2006r ) were carefully selected 
because the power spectrum of each pulsar was consistent 
with "white noise" - that is, the spectral power is statis- 
tically constant across all frequencies. However, the timing 
residuals of most pulsars (including some in our sample) are 
not consistent with white noise; rather they are affected by 
a variety of phenomena in cluding changes in th e interstellar 
medi um ||You et al.l 120071 ). calibration errors l|van StratenI 
120061 ) and irregular spind own behaviour known as "timing 

noise" (JHobbs et al.ll2010l). 

Previous authors (jLommen fc Backeill200ll : |jenet et al.l 

[2OOJ) have addressed the issue of characterising the GW sig- 
nals in pulsar timing residuals expected from SMBHBs in 
the early phase of coalescing. However, both works consid- 
ered very specific applications to known astrophysical sys- 
tems, namely the GWs being emitted from the radio galaxy 
3C66B or the Galact ic Centre (Sagittarius A* ) and nearby 
meissive dark objects. iLommen fc Backed (|200ll ) showed that 
the majcimum possible induced timing residual caused by a 
binary black hole in S agittarius A * is ar ound 14 ns, which 
is below current limits. Ijenet et al.l ([2004) showed that they 
could r ule out a propose d SMBHB system at the core of 
3C66B (|Sudou et al.ll200lj ) wi th 95% confidence using a pub- 
licly available pulsar dataset IjKaspi et al.ll 19941 ). 

The main aim of this paper is to calculate the sensitiv- 
ity of pulsar timing data to individual sources of sinusoidal 
GWs. The analysis takes account of all the issues affect- 
ing these data, including fitting of pulsar parameters, small 
amounts of non- white noise and sampling effects. This paper 
is organised as follows: § [5] describes the observations used 
to produce the sensitivity curves, §|3] describes our method 
for detecting significant sinusoids in pulsar timing residuals, 
§|3]gives our results and describes some implications and §[S] 
concludes the paper. The Appendix contains extra details of 
the detection technique we have developed. This includes a 
more thorough description of the issues encountered, in par- 
ticular those caused by fitting for pulsar parameters and the 
irregular sampling of the timing residuals. 



2 OBSERVATIONS 

The millisecond pulsar timing process usually consists of 
using a large-aperture telescope to observe a particular pul- 
sar and forming a mean pulse profile using an ephemeris 
to fold the incoming data at the correct apparent period. 
Since, on average, millisecond pulsar profiles are largely in- 
variant, a shift between the standard template pulse profile 
and the observed profile can be established, leading to a site 
arrival time. After transforming to the arrival time at the 
Solar System barycentre we obtain a barycentric TOA. Af- 
ter clock corrections are applied, and a pulsar model is fitted 



we obtain a timing residual. Timing residuals are generally 
dominated by noise, but may contain systematic errors in- 
duced by our instrumentation and more subtle effects such 
as those induced by GWs. 

The obser v ations used in this analysis were p ublished by 
IVerbiest et al] (|2008l ) and IVerbiest et all l|2009l) . who pre- 
sented results from observations of 20 pulsars using the 
Parkes radio telescopqj. Many of these pulsars exhibit some 
low-frequency noise in their timing residuals, which must be 
accounted for in our analysis. For two of the pulsars, PSRs 
J1824-2452 and J1939-I-2134, the residuals are dominated 
by low-frequency noise that complicates the spectral analy- 
sis procedures for little gain in sensitivity. Hence we remove 
them from our sample. 

The pulsars have been timed with a weighted rms resid- 
ual of ~ 0.2 — 7 lis for a period of ~ 10 years. The specifi- 
cations of each set of timing residuals are given in Table [l] 
where, in column order, we present the pulsar name in the 
J2000 coordinate system, pulse period, dispersion measure, 
orbital period, weighted rms residual, data-span and num- 
ber of recorded TOAs. For full details of TO A estimation 
and data processing, see lVerbiest et al.l (|2009l ). The timing 
re siduals for the f ull set of 20 pulsars are shown in Figure 1 
of lVerbiest et al.l (|2009l ). 

All the observations were made in the 20 cm (1.4 GHz) 
band, except for PSR J0613— 0200 for which a timing so- 
lution was obtained in the 50 cm (685 MHz) band. Obser- 
vations between 1994 and November 2002 were made with 
either one or two 128MHz-wide bands, but these data var- 
ied greatly in quality. Observations after November 2002 
were taken with a phase-coherent dedispersion syste m, the 
Caltech-Parkes-Swinburne-Recorder-2 (CPSR2; see iBailea 
[2003), over two 64 MHz- wide observing bands centred at 
1341MHz and 1405 MHz. The typical observation length 
was 1 hour. 



3 METHOD 

3.1 GW-Induced Timing Residuals 

The timing residual induced by a GW is the integral of the 
interaction between the GW and the electromagnetic pulses 
emitted by the pulsar over the path of the pulses. For our 
analysis, we assume that the GW is emitted by a binary sys- 
tem in a circular orbit. For an equal-mass binary, the lifetime 
of an SMBHB scales as (adapted from iLommen fc Backen 
I2001I) 



r = 2.2 X 10* yr 



M 



109 Mq 



-5/3 



Pnrh 



730 daysy) 



B/3 



(1) 



where M is the total mass of the system and Porb is the 
orbital perio43. For a SMBHB with M = 10^ M© and 
Poih = 730 d (which would emit GWs with a one year pe- 
riod), the lifetime is three orders of magnitude larger than 
the typical data span of pulsar timing observations. This 



" The first eight years of TOAs for PSR J1857+0943 are obtained 
from pubUcIy available d ata collected using the Arecibo radio 
telescope and presented in lKaspi et al.l 1119941 ). 
"" Note that 2Pg w = Porb , where Pqw is the period of the emit- 
ted GWs. 
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Table 1. Eighteen of the Parkes Pulsar Timing Array pulsars and their rms timing residuals from the dataset presented in lVerbiest et al.1 
l|2009l l. 



PSRJ 



Period 

(ms) 



DM 
(cm~^pc) 



(d) 



Weighted RMS 
Residual (^s) 



Span 

(years) 



No. of 
Observations 



J0437-4715 


5.757 


2.65 


5.74 


0.20 


9.9 


2847 


J0613-0200 


3.062 


38.78 


1.20 


1.56 


8.2 


190 


J0711-6830 


5.491 


18.41 


- 


3.23 


14.2 


227 


J1022+1001 


16.453 


10.25 


7.81 


1.62 


5.1 


260 


J1024-0719 


5.162 


6.49 


- 


4.22 


12.1 


269 


J1045-4509 


7.474 


58.15 


4.08 


6.64 


14.1 


401 


J1600-3053 


3.598 


52.19 


14.34 


1.14 


6.8 


477 


J1603-7202 


14.842 


38.05 


6.31 


1.92 


12.4 


212 


J1643-1224 


4.622 


62.41 


147.02 


2.50 


14.0 


241 


J1713+0747 


4.570 


15.99 


67.83 


0.20 


14.0 


392 


J1730-2304 


8.123 


9.61 


- 


2.51 


14.0 


180 


J1732-5049 


5.313 


56.84 


5.26 


3.24 


6.8 


129 


J1744-1134 


4.075 


3.14 


- 


0.62 


13.2 


342 


J1857+0943 


5.362 


13.31 


12.33 


1.21 


22.1" 


376 


J1909-3744 


2.947 


10.39 


1.53 


0.17 


5.2 


893 


J2124-3358 


4.931 


4.62 


- 


4.03 


13.8 


416 


J2129-5721 


3.726 


31.85 


6.63 


2.19 


12.5 


179 


J2145-0750 


16.052 


9.00 


6.84 


1.82 


13.8 


377 



"" There is a gap of ~11 years between the end of the data presented in lKaspi et al.l 1119941 ) and the beginning of data collected with the 
Parkes telescope. 



means no significant chirping of the GW signal will occur 
over the duration of the observations. Therefore we can cal- 
culate analytically, using Equations (1) — (7) of Jenet et al. 
(2004), the expected GW signal in the timing residuals. The 
result is that the induced timing residuals will contain two 
sinusoidal signals which can be called the "Earth term" and 
the "pulsar term". 

However, evolution of the SMBHB is sometimes a mea- 
surable effect over the timescale of the hght travel time from 
the pulsar to Earth as, for exa mple, in the evol ution of the 
proposed SMBHB in 3C66B l|jenet et all [2004 ') which re- 
sults in two distinct periodicities in the timing residuals. 
In this work we ignore this longer timescale evolution, so 
the GW-induced quadrupolar space-time distortions at the 
Earth and the pulsar will always have the same frequency. 
However we have allowed the two periodicities to be offset 
in phase which can alter the amplitude of the signal in the 
timing residuals. We hence reduce the problem of detecting 
GW emission from a non-evolving circular binary system 
to identifying the presence of a significant sinusoid in the 
timing residuals. To confirm that a significant sinusoid in 
the timing residuals of a given pulsar is caused by GWs, one 
would need to ensure that the expected signature of the GW 
(see, e.g., iDetwcilcr 19791 ) is present in the timing residuals 
of other pulsars. 

To determine the residuals a particular SMBHB will 
induce in our data, we begin with the expe cted GW strairjj 
emitted by a single SMBHB (|Thorndl 19871 ): 



/..=4./^(^^' 



VV3 



5 C4D(2) 



[Kf{l + Z] 



i2/3 



(2) 



^ By convention h^ is the strain from a GW background while 
hs gives the GW strain from a single source. 



where Mc = [MiM2f'^ (A/i + M2) ^'^ is the chirp mass 
of the SMBHB with member masses M\ and M2, G the 
gravitational constant, c the vacuum speed of hght, / the 
observed GW frequency (which is in general different to the 
emitted frequency), z the redshift of the SMBHB and D{z) 
is the comoving distance to the SMBHB, given by 



D{z) 



c 



dz' 



(3) 



where Ho is Hubble's constant, taken to be 72 km s ^ Mpc ^ 
and E{z) = H{z)/Ho = JQ. k + ^m(l + zY under a ACDM 
cosmological model (see, e.g.. lStavridis et aU 20091). For this 
work we assume f^A = 0.7 (see, e.g., iKomatsu et al.l 120091 ) 
giving Sim — 0.3. 

A non-evolving GW source will induce a sinusoidal 
varia tion in the pulsar timing residuals, with amplitude 
(We n et al.|[2O10l ) 



A, 



-(l + cos9)sin(2<?!))sin 



ujDp(l — cos 9) 



2c 



(4) 



where uj = 27r/PGW is the GW frequency in rad s~^, 6 is the 
angle between the direction from which the GWs emanate 
and a vector from the Earth to the pulsar, </!> is the GW 
polarisation angle and Dp is the distance to the pulsar. 

An important feature of any pulsar timing analysis is 
the process of parameter determination for the model of the 
pulsar. This process is equivalent to fitting out a range of 
signals from the time series of residuals. Figure [l] shows the 
effect this can have on GW detection - a GW signal with a 
period of one year (top left panel) will be almost completely 
removed after fitting (top right panel) because this signal 
mimics an error in the pulsar position. However, a GW signal 
with a period of two years (bottom left) is only slightly at- 
tenuated by fitting (bottom right). To determine the post-fit 
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Arrival Time 



Arrival Time 



Figure 1. Attenuation of timing signals caused by pulsar parameter fitting. In each plot the abscissa is the centred MJD and the 
ordinate gives the magnitude of each timing residual; all plots have the same scale. The dotted lines indicate zero residual. The data 
plotted are formed by adding a simulated GW signal to the timing residuals for PSR J1909— 3744 which are described in Section [2] 
and in IVerbiest et al.l 1 120091) . The top row shows a GW signal with a period of one year (top left) being completely removed after 
fitting (top right). The bottom row shows a GW signal with a period of two years (bottom left) being largely unaffected by the fittirig 
procedure (bottom right). These figures were produced using simulated data from the pulsar timing package Tempo2 tHobbs et al.ii200a 
lEdwards et alll2006l: iHobbs et aLll2009h . 



timing residuals, we add tlie effect of a sinusoid al GW point 
source dire ctly to the TOAs using Tempo2 (JHobbs et al.l 
120061 . 120091 ) and then perform the standard pulsar timing 
fitting procedure on these modified TOAs. 



3.2 Producing the Sensitivity Curve 

The detection of a sine wave in the presence of noise with 
known statistics is a well-studied problem with a simple opti- 
mal solution, the maximum likelihood estimator. A number 
of algorithms can be used, depending on the characteris- 
tics of the data. In this case the problem is complicated by 
the fact that the data are irregularly sampled and the noise 
consists of at least two components. The noise has a white 
component which varies from sample to sample; this compo- 
nent is well understood and we have a variance estimate for 
the white noise at each sample point. The noise also has a 
non- white component for which the source is unknown. We 
assume that it has a smoothly varying power spectrum and 
attempt to estimate this from the data. 

We use one of the most common spectral esti- 
mation tools: an unweighted Lomb-Scargle periodogram 
IjPress et al.lll992l ). By "unweighted", we mean that the in- 
dividual TOA errors are not taken into account when calcu- 



lating the power spectrum. The Lomb-Scargle periodogram 
technique is not valid for datasets which exhibit a steeply 
sloping spectrum (the timing residuals of some young pul- 
sars exhibit very high power levels only at low frequencies). 
While many of our timing residuals do not exhibit a flat 
Lomb-Scargle periodogram, none of our 18 datasets has a 
sufficiently steeply sloping spec trum to invalidate th is ap- 
proach. It has been argued (see lCumming et al.lll999r ) that 
a 'fioating mean' periodogram should be used to obtain the 
correct spectral estimates when detecting long period sig- 
nals in sparsely and unevenly sampled data; this potential 
improvement will be addressed in a future paper. We briefly 
describe our approach for producing a sensitivity curve here; 
full details are provided in the Appendix. 

To make a detection of a significant sine wave in our 
timing residuals, we make a simple model of the noise in 
the Lomb-Scargle periodogram of the timing residuals and 
use this model to define a set of detection thresholds. These 
thresholds are set high enough that the probability of record- 
ing a detection at any frequency across the entire observed 
power spectrum when no signal is present is 1% (the "false 
alarm probability"). 

We then inject simulated GW signals of known fre- 
quency with random polarisations, random sky locations 
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and at a range of strain amplitudes to determine the strain 
that gives a detection in the corresponding frequency chan- 
nel in 95% of the simulations. The GW-induced quadrupo- 
lar space-time distortions at the Earth and the pulsar can 
interfere constructively or destructively. This process gives 
the sky- and polarisation-averaged sensitivity as a func- 
tion of GW frequency over the range / ~ (10 years) ^^ to 
/ ~ (1 month)" . 

There are two aspects to this detection strategy, namely 
the false alarm probability (1%) and the probability of mak- 
ing a detection (95%). Using a false alarm probability of 1% 
means that in a particular simulated dataset, any detection 
made will be a 3-a detection. However, the sensitivity of a 
pulsar GW detector is not the same for every possible GW 
source; for instance, a single pulsar cannot be used to detect 
GWs propagating along the line of sight from the Earth to 
the pulsar. Our sensitivity curves give the GW amplitude at 
which the probability of making a 3-a detection at a ran- 
dom position on the sky and with a random polarisation 
is 95%. If the GW polarisation and the pulsar-Earth-source 
angle were favourable (e.g. d = n/2 and (j) = 7i"/4 in equation 
[4]), then for a single pulsar an improvement by a factor of 
~ 10 — 15 in sensitivity could be achieved compared with the 
sky-averaged sensitivity. One of the advantages of timing an 
array of pulsars to detect sinusoidal GWs is that there is a 
significantly smaller region of sky in which a full array has 
low sensitivity. 

We are interested in answering two questions. The first 
is "What is the largest GW source at a particular frequency 
that could be present in the measured timing residuals?" 
This will give an upper bound on the amplitude of individ- 
ual GW sources in our data at that frequency. This question 
is answered by comparing simulated GW sources to our ob- 
served timing residuals. We simulate GW sources at a given 
frequency at random sky locations and adjust the amplitude 
of these sources until the power of the GW sinusoid exceeds 
the power in the observed timing residuals at that frequency 
in 95% of simulations. This approach gives the most conser- 
vative upper limit, since we are allowing for the possibility 
that all the power we observe at this frequency results from 
one sinusoidal GW point source. We will determine this up- 
per limit for our datasets in § 4. 

The second question is "If there were a GW source with 
a particular frequency somewhere on the sky, what is the 
minimum strain amplitude that would produce a detectable 
signal at that frequency in our dataset?" In order to answer 
this we add simulated sinusoidal GW signals to our TOAs, 
perform the standard pulsar timing analysis and calculate 
the minimum amplitude at which we would detect a sig- 
nificant sinusoid at the input GW frequency in our data if 
we had collected that dataset at a telescope. This means we 
must account for all the sources of noise in our pulsar de- 
tector. The threshold for detection at any frequency across 
the observed power spectrum will often be ~3 times greater 
than the locally-averaged power level. This gives our sensi- 
tivity to detecting these sinusoids, rather than just limiting 
their amplitude. For very large amplitude sine waves with 
period > Tobs, a signal will often be detectable at a slightly 
higher frequency than the input frequency because we can 
detect the side lobes of the large input signal. We have not 
allowed detections at different frequencies to the input GW 



frequency in our implementation. The sensitivity curve cor- 
responding to our datasets is derived in § 4. 



The periodogram frequency range is from 



to 



Tabs 2T„b, 

for a single pulsar (where Tobs is the time-span of the obser- 
vations and A'^pts is the number of timing residuals for that 



pulsar). Note that 



would be the Nyquist frequency for 



that pulsar if its timing residuals were regularly sampled. If 
we have multiple pulsars then we can perform a weighted 
sum of their power spectra to increase our sensitivity. For 
this work we first make a simple frequency-dependent model 
of the noise in the pulsar power spectrum and then weight 
each pulsar by the inverse of the noise model for that pulsar. 
For the case of spectrally white timing residuals, the noise 
model will be independent of frequency, so this is equivalent 
to weighting by the inverse variance of the timing resid- 
uals. To perform the sum requires the use of a common 
frequency gridding, so when analysing multiple pulsars the 
periodogram ranges from (30 years) ~^ to (4 weeks) ~^ for all 
pulsars. 

Our detection technique is straightforward to imple- 
ment, but for many pulsars with differing data-spans and 
noise properties etc., it will not be optimal. After adding 
a simulated GW signal to each dataset, some of the power 
at the frequency of the GW signal will be leaked into adja- 
cent channels, meaning that the noise model will be higher 
near the GW frequency, leading to fewer detections. We also 
use a simple weighting scheme to combine multiple pulsars 
which gives a factor of ~ 5 improvement over a simple, non- 
weighted addition of the power spectra for the different pul- 
sars. However, the exact weighting used in an incoherent 
detection scheme such as ours does not significantly change 
the overall sensitivity of the array. A small improvement in 
sensitivity may also be gained by allowing for the evolution 
of the GW source over the light-travel-time from the pulsar 
to the Earth and then searching for a two- frequency response 
in each pulsar's power spectrum. The 18-pulsar array sen- 
sitivity could also be improved by 'phasing up' the timing 
array to enable a coherent sum of the GW signal in each 
dataset. However, such a detection scheme is considerably 
more complex and will be addressed in a future paper. 



4 RESULTS AND DISCUSSION 

In this section we present the sensitivity of the Parkes Pulsar 
Timing Array (PPTA) to sinusoidal GW sources using the 
dataset described in Table[T]and account for all the observed 
features in the sensitivity curves. We also describe some of 
the implications of the non-detection of sinusoidal GWs in 
our dataset and give predictions for a future timing array 
project using the SKA. 



4.1 The sensitivity using some individual pulsars 

In Figure [2] we plot the sky- and polarisation-averaged 
sensitivity curves for PSRs J0437— 4715 (thin solid line), 
J1713-H0747 (dashed line) and J1857-H0943 (dot-dashed 
line) where each pulsar has been analysed individually. 
The open triangles on the plot indicate that the plotted 
"detectable" amplitude at that frequency value is a lower 
bound. The thin dotted line indicates the sensitivity of PSR 
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Figure 2. Sensitivity curves for PSRs J0437-4715 (thin solid line), J1713+0747 (dashed), J1857+0943 (dot-dashed) and the 18-puIsar 
timing array assuming an incoherent detection scheme is used (thick solid line) . The abscissa gives the GW frequency, the ordinate gives 
the minimum detectable strain amplitude of a sinusoidal GW point source with a random polarisation, phase and sky-position. The 
thin dotted line is the maximum attainable sensitivity using PSR J0437— 4715 assuming optimal sky-location and polarisation of the 
GW source. An open triangle indicates that the plotted value is in fact a lower bound on the detectable amplitude at that frequency. 
The straight triple-dot-dashed lines indicate the expected signal from an individual SMBHB with equal member masses of 10 Mq or 
10^*^ Mq if it were located at the mean distance of the Virgo cluster. The 'x' symbols are the expected signals at the Earth now and 
at PSR J1857+0943 ~2700 years ago caused by the proposed SMBHB at the core of the radio galaxy 3C66B. The '*' symbol is the 
expected signal caused by the proposed SMBHB a t the core of OJ287 . The '-I-' symbol is the GW strain and frequency emitted by a 
typical resolvable SMBHB as plotted in Figure 2 of lSesana et al.l 1 120091) . Also shown on the plot is the upper limit on the amplitude of 
sinusoidal GW point sources as a function of frequency using the 18-pulsar timing array (thick dotted line); in this case the ordinate 
gives the maximum amplitude GW source which could be present in our data. 



J0437-4715 to a hypothetical SMBHB located at a right as- 
cension of 4*^37™ and a declination of -(-42°45™ and emitting 
purely 'plus' polarised GWs. This line indicates the max- 
imum sensitivity obtainable with this dataset if the GW 
source position and polarisation were favourable in every 
simulation. The ratio of this thin dotted line to the thin 
solid line gives the factor of ~ 10 — 15 improvement in sen- 
sitivity for optimal sky-location and polarisation discussed 
in § 13.21 Also shown are the expected signals at a range of 
frequencies from two hypothetical SMBHB systems at the 
mea n distance of th e Virgo cluster (taken to be 16.5 Mpc, 
from lMei et al.ll2007h . with equal member masses of 10^ M0 
or IO^^Mq. 

The reduction in sensitivity caused by fitting for the 
pulsar's position will be at the same frequency of (lyr)~^ 
for all pulsars. However, fits for orbital parameters will also 



reduce sensitivity to GWs, but at difi'erent frequencies for 
each pulsar. All pulsars exhibit a reduction in sensitivity at 
low frequencies; this is caused by the fit of a quadratic poly- 
nomial to the TO As required to model the pulsar spin-down, 
as well as the fitting of "jumps" to many of the datasets to 
connect the timing residuals obtained with different backend 
systems (see below). 

As the GW frequency increases, the strength of the sig- 
nal in our residuals becomes weaker for a given strain, as 
described by equation (|3|. At the highest frequencies, our 
sensitivity is limited by the sampling of the timing residuals. 
This is particularly evident in the sensitivity curve for the 
18-pulsar timing array where there is a turn-up in the sen- 
sitivity curve at the last few frequency values corresponding 
to a decrease in sensitivity there. 

The sensitivity of our detection technique to low- 
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frequency sinusoidal GWs in irregularly-sampled data 
(where the GW period is similar to the data-span) is re- 
duced compared to treating regularly-sampled data since 
there is no clear way to distinguish between the excess low- 
frequency noise seen in many millisecond pulsars and spec- 
tral leakage from the low- frequency GWs. Some pulsars in 
our sample do not exhibit excess low-frequency noise (e.g., 
PSR J1857+0943), meaning that the power spectrum with 
no GWs added may be modelled with a constant. However, 
as soon as a low-frequency sinusoidal GW source is added 
to these residuals, leakage from the low-frequency signal is 
difficult to distinguish in our technique from standard pul- 
sar timing noise and interstellar medium variations, so our 
model of the power spectrum must account for this confu- 
sion. In a regularly-sampled time series with weak red noise, 
spectral leakage is less severe and thus there is no such con- 
fusion. 

In the sensitivity curve for PSR J0437— 4715 there is a 
loss of sensitivity at a frequency of (540 days) ~^, or ~21 uHz. 
This is caused by the fitting of several constant time offsets 
between the data collected using different observing backend 
systems; such offset fits absorb GW power at low frequen- 
cies. If overlapping data exist between the different observ- 
ing backends, these offsets can be precisely determined and 
held fixed in subsequent processing. Even if no overlapping 
data exist, it is sometimes possible to eliminate these arbi- 
trary offsets without losing phase connection in the timing 
solution. Our analysis takes i nto account all of the offsets 
fitted bv lVerbiest et al.l (120091 ). There is also a loss in sensi- 
tivity just above the (lyr)~^ frequency for this pulsar. This 
is caused by the sampling of the dataset. 



4.2 The sensitivity of the PPTA and some likely 
single sources 

Figure [5] contains the sky-averaged sensitivity attainable 
for the 18 pulsars in our dataset assuming an incoher- 
ent detection scheme is used and the GW source position 
and polarisation are unknown. The plotted frequency range 
(30 years) "'^ - (4 weeks) ~^ is chosen to demonstrate the 
high- and low-frequency sensitivity limits for our pulsar tim- 
ing datasets. At the lowest frequencies, our sensitivity is 
limited by the fact that our longest dataset is much shorter 
than 30 years and by the necessary period derivative and 
jump fits. At the highest frequencies we are limited by the 
sampling of our timing residuals; that is, (4weeks)~^, the 
nominal Nyquist frequency for the PPTA. 

Figure [2] also show s the upper limit att ainable using 
the 18 pulsars from the IVerbiest et al.l (|2009|) dataset (the 
thick dotted line at the bottom). This limit curve was ob- 
tained with 95% confidence as described in Section [3.21 and 
in the Appendix. For some pulsars a different-order poly- 
nomial model to the detection case was chosen in order to 
accurately model the pow er spectrum with no GWs added. 
iLommen fc Backeii (|200ll ) placed a 99% confidence limit 
showing that they could rule out signal amplitudes as small 
as 150 ns in their residuals at a period of 53 days, correspond- 
ing to SMBHB orbital periods of 106 days. Using our longer 
datasets and the same 99% confidence level, we can place 
a better limit of around 120 ns at this frequency. At signal 
periods of 1000 days where some of our datasets exhibit ex- 
cess low-frequency noise, we obtain a 99% confidence limit 



of 190 ns which is worse than the lLommen fc Backeil (|2001h 
limit of 170 ns. However, there is no evidence that their anal- 
ysis takes into account the effects of red noise present in their 
residuals. 

The two 'x' symbols in Figure [2] indicate the expected 
strain amplitude and frequency of th e proposed SMBH B at 
the core of the radio galaxy 3C66B (jSudou et al.l 120031 ). In 
order to determine the expected strain amplitude, we use 
equation ([2]) with the redshift and masses given in the orig- 
inal paper (mi = 4.91 x lO^"M0,m2 = 4.91 x 10*^1^10,2 = 
0.0215) and a distance to the source of 90 Mpc, implied 
by the low-redshift distance approximation D = cz/Hq. 
The frequencies of the signal at the Earth and at PSR 
J1857-h0943 (/Earth =_2Za88yr, ^1857-^0943 = l/6.24yr) 
were obtained from Ijenet et al.l I 2004h. This system was 
ruled out with 95% confidence by Ijenet et al.l (|200J). Our 
results show that even with a blind search of the Verbiest 
et al. data, where we know neither the sky position nor the 
frequency of the GWs, we would detect the GW-induced 
oscillations at the Earth caused by this source. The ex- 
pected signal is well below the plotte d sensitivity curve for 
PSR J1857-f0943 even though i Jenet et al., (2004) only used 
the publicly available timing residuals for PSR J1857+0943. 
However, their technique is analogous to our limit technique, 
whereas the sensitivity curve plotted for PSR J1857-)-0943 
in Figure [2] assumes we are aiming to detect such sources 
of GWs. Furthermore, our sensitivity curve is sky-averaged 
whereas they used the known position and frequency of the 
proposed GW source in their analysis (by chance it had a 
very favourable sky-location with an angle of 81.5° between 
the Earth- pulsar vector and the Earth-3C66B vector) [j 

The '*' symbol in Figure [2] indicates the expected GW 
strain and frequency for the candidate SMBHB in the blazar 
OJ287. A ~12yr- periodic signal has been identified in its 
optical outbursts dSillanpaa et al.l 1996), but other parame- 
ters of the system are not well-constrained. We parametrise 
the SMBHB as follows: member masses 1.3 x 10* M© and 
1.8 x 10^" Mq, orbital period 9 years (observed GW period 
4.5 years), eccentricity Qj, redshift 0.306, distance 1.3 Gpc. 
The distance was again obtained using D = cz/Ho, which 
is an acceptable approximation given the imprecision in the 
other parameter measurements and the fairly low redshif t 
of this system (see footnote 1 in lDavis fc Lineweaverll2004l ). 
The GW signals emitted by this system induce timing resid- 
uals of around 5 ns which a re well below current limits. 

In ISesana et al.l l|2008i ) a study was presented of the 
generation of the stochastic gravitational-wave background 
from the cosmic population of SMBHBs. This work showed 
that the stochastic background of GWs is likely to be de- 
tected using a pul sar timing array in the near future. In 
ISesana et al.l ((20091) the individual resolvable SMBHBs were 
considered. They predicted that at least one SMBHB will in- 
duce timing residuals around 5 — 50 ns, which is below our 
current sensitivity. We choose (from the upper left panel 
of their Figure 2) a representative resolvable single source 
from their simulations, with an emitted GW frequency of 



^ Ijenet et al.l 112004) also underestimated the distance to the pro- 
pose d GW source in 3C66 B by around 10%. 

* In lValtonen et al.l 1 120091) the eccentricity is estimated to be 0.7, 
but we do not consider eccentric SMBHBs in this paper. 



8 Yardley et al. 




Frequency (Hz) 

Figure 3. Sensitivity of the PPTA using the 18-pulsar Ver- 
biest et al. dataset for detecting signals from SMBHBs located 
at the sky-position and mean distance of the Virgo cluster. The 
abscissa gives the GW frequency, the ordinate gives the minimum 
detectable strain amplitude of a sinusoidal GW point source em- 
anating from the direction of the mean sky-position of the Virgo 
cluster with a random polarisation and phase. The open triangles 
indicate that the plotted value is a lower bound on the detectable 
amplitude at those frequencies. The dot-dashed lines indicate the 
expected signals from three different types of SMBHB if they were 
located in the Virgo cluster, with equal member masses 10^ Mq, 
lOl°M0 and lO^M© as labelled. 



2 X 10~* Hz and a characteristic induced timing residual of 
25 ns. The signal from this source is indicated by the '-I-' 
symbol in Figure [2] This is a typical resolvable SMBHB, 
thus it is likely that several sources will emit GWs with a 
larger amplitude than this. We emphasise that we do not 
yet have long data-spans with sufficiently low rms residual 
to detect such sources. A stochastic background of GWs may 
be detected using a pulsar timing array within the next few 
years. 

The formation of SMBHBs is more likely in galaxy clus- 
ters. The nearest galaxy cluster to Earth is the Virgo cluster. 
In Figure [3] we examine the possibilities for pulsar timing to 
detect GWs generated by SMBHBs in the Virgo cluster. The 
mean sky-position of this cluster i s at a right asc ension of 
12'' 30" and a decfination of 4-12° (|Mei et al.ll2007l ): to pro- 
duce this sensitivity curve all simulated GW signals come 
from this direction. The plotted sensitivity curve indicates 
that, with a false alarm probability of 1%, we have a bet- 
ter than 95% probability of detecting sinusoidal signals in 
our timing residuals caused by 10^° M© - 10^° M© SMB- 
HBs in the Virgo cluster with any polarisation at a range 
of frequencies and marginally also some 10^ M© — 10^ M© 
SMBHBs. 

The PPTA sensitivity is complementary in GW fre- 
quency to the LIGO, VIRGO and LISA sensitivities. In Fig- 
ure |4] we give the detection sensitivity of some current and 
future GW detection experiments. Also shown on the plot 
are some likely sources in each of the detectable bands. This 
sensitivity curve now almost covers the full GW frequency 
range from ~nHz through to ~mHz; this frequency cover- 
age will enable the study of the evolution of GW-emitting 
systems. 

To obtain the LISA sensitivity curve, we have assumed 



the standard parameters for the LISA design and that it 
aims to detect sources at a signal-to-noise ratio of three. 
The LIGO sensitivity curves are obtained from the stated 
design goals of the project. 



4.3 The implied constraint on the merger rate of 
SMBHBs 

Non-detection of single-source GWs in the Verbiest et al. 
data enables an upper limit to be placed on the rate of 
super- massive black hole me rgers 11 Wen et al.l 2010l ). To use 
the techniques presented in IWen et al.l (2010|), it is neces- 
sary to calculate the limiting sensitivity of our array at a 
matrix of GW frequency and strain values. The frequency 
values chosen were 50 logarithmically-spaced frequencies be- 
tween (30 years) ^ and (4 weeks)" , while the strain values 
were 50 logarithmically-spaced amplitudes between 10"^'' 
and lO"^''. The sensitivity at a frequency of (lyr)~^ was 
also calculated, resulting in 51 frequency values overall. 1000 
Monte Carlo iterations were used at each value of GW fre- 
quency and strain. The sensitivity matrix obtained gives a 
95% confidence contour which is consistent with the 95% 
confidence upper bound obtained earlier (the thick dotted 
line in Figure [5}. 

This sensitivity matrix is used to provide an upper limit 
on the differential rate of SMBHB coalescence per logarith- 
mic redshift and chirp mass. The results of this analysis 
are shown in Figure [S] Our data do not yet con s train the 
merging framework s discussed by I Jaffe fc Backer! (12003) or 
(jSesana et al.l 120081 ) at the range of chirp masses we have 
considered. However, in coming years some of the high-mass 
and high-redshift predictions may be ruled out or confirmed 
using pulsar timing. 



4.4 A predicted SKA sensitivity curve 

Figure 2] also gives a predicted sensitivity curve for the 
Square Kilometre Array (SKAJj. To produce this figure we 
chose 100 pulsars from the Australia Telescope National Fa- 
cility pulsar catalogue IjManchester et al.ll2005^ ■ We have as- 
sumed we can time each pulsar with an accuracy of 20 ns 
over five years, obtaining one timing point per pulsar every 
two weeks. We have also assumed that their power spectra 
will be statistically white. This means the plotted sensitiv- 
ity is a lower bound on what is achievable with the SKA 
for the assumed parameters, especially at low frequencies 
where we expect higher noise levels caused by the stochastic 
background of GWs and intrinsic pulsar timing noise. 

The simulated SKA data are regularly sampled with 
equal error bars, which means that the level of spectral 
leakage will be much lower than that observed in irregu- 
larly sampled datasets with highly variable error bars. This 
means the confusion between red noise and low-frequency 
signal is no longer an issue in these simulations because a 
sinusoidal GW signal will induce a very narrow peak in each 
pulsar's power spectrum, even at low frequencies. We have 
therefore modelled each pulsar power spectrum with a con- 
stant. 



See |http://www.skatelescope.org7] 
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Figure 4. Sensitivity of some current and future GW observatories to individual GW sources as a function of frequency. The abscissa 
gives the GW frequency, the ordinate gives the minimum detectable strain amplitude of a sinusoidal GW point source with a random 
polarisation, phase and sky-position. For the pulsar timing array sensitivity we assume an incoherent detection scheme is used and the 
GW source position is unknown. The open triangles indicate that the plotted sensitivity at that frequency is a lower bound. The plot 
also shows some potentially detectable sources in the three frequency bands. The straight lines indicate the expected signals from two 
different types of SMBHB if they were located in the Virgo cluster, with equal member masses lO^M© and lO^'^Mo as labelled. The 
'x' symbol is the expected signal at the Earth caused by the proposed SMBHB at the core of the radio galaxy 3C66B. The '*' symbol 
is the expected signal caused by the candidate SMBHB a t the core of OJ287 . The '+' symbol is the GW strain and frequency emitted 
by a typical resolvable SMBHB as plotted in Figure 2 of ISesana et al.l 1120091 ). "Unresolved galactic binaries" include white-dwarf and 
neutron-star binaries. "Coalescing binary black holes" show the expected range of signals from the final inspiral of black-hole binary 
systems. The "Current" LIGO sensitivity shows the capabilities of existing datasets, while "Advanced" LIGO expects to improve GW 
sensitivity by two orders of magnitude. "SN [supernova] core collapse" and "NS-NS [neutron star] coalescence" are typical signals that 
LIGO expects to detect. 



There are three prominent losses in sensitivity - at fre- 
quencies smaller than (Tobs)~ and at periods of one year 
and six months. The partial loss in sensitivity at a period of 
six months (~ 6 x 10~* Hz) is caused by fitting for the pulsar 
parallax. The total loss in sensitivity at GW periods of one 
year could be mitigated using independent measurements 
of the position of the pulsar, for example using very-high- 
precision interferometry; such precision may be available in 
the SKA era. 

The SKA sensitivity curve shown in Figure |4] is calcu- 
lated assuming we do not know the location or frequency of 
a potential GW source; using these two additional pieces of 
information it will be possible to confirm or deny the bina- 



rity of the massive dark object at the core of OJ287, as well 
as resolve many of the SMBHBs predicted by ISesana et al.l 
(2009). Using the SKA and LISA, it will also be possible to 
observe the full evolution of some SMBHBs from emitting 
GWs in the pulsar timing band (during the early phases 
of coalescenc e) to emitting GW s in the LISA band (during 
coalescence) jPitkin et al.ll200a ). 



5 CONCLUSION 

We have presented the strain sensitivity of the Parkes Pul- 
sar Timing Array to sinusoidal point sources of GWs as a 
function of frequency. The sources most likely to produce a 
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Figure 5. Upper limit on the rate of SMBHB mergers as a function of redshift (witti a range of chirp masses) given the non-detection 
of any GW sources in the Verbiest et al. datasets. The open triangles give the upper limit on the SMBHB merger rate for the Verbiest 
et al. dataset and t he open squares give t he limit for the simulated SKA datasets. The shaded region indicates the expected coalescence 
rate obtained from ljaffe &: Backer! 1 120031 ) as well as data from the Sloan Digital Sky S urvey llWen et al.ll2009l ) for SMBHB systems of 
chirp mass as labelled in each panel. The dashed line indicates the merger rate based on lSesana et al.l 1 120081) . 



detectable sinusoid in the pulsar timing frequency range are 
super-massive black-hole binary systems in the early phases 
of coalescence at the cores of merged galaxies. The sensi- 
tivity curve is analogous to the LIGO, VIRGO and LISA 
sensitivity curves and indicates the unique GW frequency 
range accessible with pulsar timing. These results can be 
used to place an upper bound on the number of coalescing 
binary systems of a given chirp mass as a function of red- 
shift. Current observations do not yet rule out any likely 
GW sources. 
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APPENDIX: DETAILS OF THE DETECTION 
TECHNIQUE 

In this section we give a detailed description of our detection 
technique, in particular describing some of the problems that 
arose during this treatment. 

Our Technique for Producing a Sensitivity Curve 

Our method for creating curves showing the sensitivity of 
our timing residuals to GW-induced sinusoidal signals from 
individual SMBHBs takes into account the non-Gaussian 
noise which is a feature of many timing residual datasets. 
To produce a sensitivity curve for a given set of pulsars and 
their timing residuals, we use a 3-step process as follows: 

(i) We choose logarithmically spaced GW frequencies be- 
tween 7f^ and 7^ 



(single pulsar) or between (30 years) ^ 
and (4 weeks)" (multiple pulsars). The frequency sampling 
we used for multiple pulsars requires over-sampling each 
power spectrum by a factor 30yr/Tobs for that pulsar, 
(ii) At each frequency, we: 

(a) add the effect of a sinusoidal GW point source with 
angular frequency 27r/i, amplitude hs and random sky- 
position and polarisation to the TOAs, as described in 
equation Q. 

(b) process the data using the Tempo2 pulsar timing 
software to obtain post-fit timing residuals. 

(c) run a detection algorithm (described below) on the 
post-fit residuals which reports either a detection or a 
non-detection. 

a large number of times 



(d) repeat steps (ii)(a) - (ii)(c 



(we use 1 000 iterations) and record the detection percent- 
age. 

(e) If we have detected (95 ± 1) % of the signals then we 
have satisfied our detection criterion and we record fi and 
hs, which places a point on the pulsar timing sensitivity 
curve. If this criterion is not satisfied, adjust hs higher if 
too few detections have been made and lower if too many, 
then return to step |(ii)(a)[ 

(iii) Select the next frequency in the grid and repeat. 

Our detector functions as follows: 

(i) For each pulsar in the input data, we calculate a 
power spectrum of the residuals using a Lomb-Scargle peri- 
odogram, with the frequency range described above. 

(ii) We smooth the power spectrum by taking the loga- 
rithm of the power values and using a boxcar median filter 
(by default the number of points in the filter is 11 times 
the oversampling factor for that pulsar, in order to account 
for the correlated spectral estimates induced both by over- 
sampling and by the irregular time sampling of the timing 
residuals). 
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(iii) We use a least-squares fit to obtain a low-order poly- 
nomial (i.e. of order less than six) that provides a simple 
model of the median-smoothed log-spectrum. The median- 
smoothing and model-fitting are performed only on those 
points for which the frequency is ^ (T'obs)^ • This three- 
step spectral modelling process ensures that the simulated 
GW source is not included in the model as part of the 
noise in the spectrum; this is particularly important at the 
low- and high-frequency edges of the power spectrum. When 
analysing the data collected from multiple pulsars we com- 
bine their power spectra using a weighted sum. The weight 
used for each pulsar is the inverse of the simple frequency- 
dependent model of the power spectrum for that pulsar. 

(iv) We multiply the noise model obtained above by a fac- 
tor of ~2-3 determined from simulation (see below) to define 
a set of detection thresholds for any given false alarm prob- 
ability (we use Pf — 1%). These detection thresholds are 
set such that the probability of any observed power across 
the whole spectrum being greater than the threshold when 
there is no signal present is 1%. 

(v) If the measured power in the channel containing the 
input GW frequency is greater than the detection thresh- 
old in that channel, then we have made a detection of a 
significant sinusoid. 

Some of the simulated sinusoidal GW point sources pro- 
duce large signals in the timing residuals, depending on their 
polarisation and location on the sky. If a set of timing resid- 
uals showed evidence of a strong signal, a typical analysis 
would use a model of the pulsar with the fewest possible 
parameters (i.e. a period, period-derivative and any arbi- 
trary phase offsets) to obtain residuals and then examine the 
dataset more closely. To simulate this process, in step (ii)(b) 
above, we calculate the full parameter fit as normal, but if 
the reduced-x^ is larger than 20, then we instead only fit for 
the pulsar period, spin-down and jumps between datasets. 



Our Technique for Producing an Upper Limit 



As described in Section 13.21 our technique for ruling out 
GWs with a particular strain amplitude as a function of 
frequency is much more straightforward than attempting to 
make a detection of the sinusoid induced by GWs emanating 
from SMBHBs. The important assumption in producing a 
correct upper limit without assuming anything about the 
statistics of the data in question is that, at any frequency 
in our power spectrum, the power caused by GWs cannot 
be more than the observed power; otherwise, we would have 
observed a higher power level at that frequency. That is, we 
assume that all the power at a given frequency is caused by 
GWs and then calculate the GW strain which gives a power 
greater than this level 95% of the time. 

To produce this limit, we first calculate the power spec- 
trum of the observed timing residuals of each pulsar us- 
ing the Lomb-Scargle periodogram. We then make a sim- 
ple polynomial model of the noise in this spectrum and use 
the inverse of this noise model as the weight in calculat- 
ing a weighted sum of the power spectra. This weighted 
and summed spectrum is the detection threshold. We then 
simulate noiseless GW signals (which are manifested as a 
pure sinusoid in each dataset), fit out as many of the pul- 
sar parameters as possible from each sinusoid and calculate 



the same weighted sum described above (that is, using the 
noise model calculated for the observed timing residuals). 
Comparing this weighted sum of sinusoids to the detection 
threshold, we can scale the strain amplitude until we can 
detect the signal in 95% of detection attempts. We can then 
rule out the existence of any stronger GW sources at this 
frequency (with random sky position and polarisation) with 
95% confidence. 

The False Alarm Probability 

We used simulation to calculate the correct detection thresh- 
old for a given dataset and a false alarm probability of 1% 
across the whole spectrum. Any detection made will thus 
be a 3-cr detection. The statistics of each channel in the 
power spectrum approximately follow a x^-distribution, but 
many other issues change the statistics of each channel, as 
described below. 

As soon as we add a large GW signal to our data in 
channel i, the statistics of channel i follow a non-central x^- 
distribution or a Ricean distribution. This does not affect 
the false-alarm probablility determination but would affect 
analytical determinations of pulsar timing sensitivity. 

Other effects that change the statistics of each spec- 
tral channel include the irregular sampling of the time series 
(which can cause correlated estimates of the power at some 
frequencies), the oversampling of the power spectrum when 
analysing multiple pulsars (which means that the peaks in 
the power spectrum will be more fully resolved and thus the 
peak value is higher) and the median filtering (which lowers 
the height of each peak in the spectrum as well as raising 
the troughs). 

Our method for calculating the height of the 3-cr detec- 
tion threshold was to simulate many realisations of white 
noise with an rms of 100 ns and the same sampling as the 
original time series. Then, without performing any of the 
pulsar parameter fits or adding the effect of a SMBHB, 
calculate the average detection rate for any peak in the 
power spectrum to be greater than some estimated detection 
threshold. To perform the "detection" in this case we sim- 
ply find the mean of the power spectrum (since the timing 
residuals are consistent with white noise) and then make the 
estimated threshold a factor of '-^ 2 higher than this mean. 
This factor is adjusted until the average detection rate for 
detections being made anywhere in the observed power spec- 
trum equals the false alarm probability. In general the de- 
tection threshold had to be set at a factor of 1.3—2.5 higher 
than the threshold implied by assuming that each spectral 
channel follows a x^-distribution. 

Modelling the Power Spectrum 

In Figure [6] we show a sample dataset with a very low fre- 
quency GW source injected and the models used for the 
three individual pulsars whose sensitivity is displayed in Fig- 
ure [2l In general, the models chosen are conservative in the 
presence of red noise to minimise the number of spurious de- 
tections at low frequencies. These figures demonstrate some 
general features of the power spectral models used. In par- 
ticular, the models account for the varying levels of red noise 
and the possibility of signal leakage. When limiting the am- 
plitude of the single sources that could be present in our 
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Figure 6. A typical power spectrum of each of three timing residual datasets, where we have added a large amplitude low-frequency 
source to each set of timing residuals. The abscissa gives the frequency examined, the ordinate gives the power in arbitrary units including 
constant offsets applied to the spectra of PSRs J1857+0943 and J1713-I-0747 to separate the spectra in making this plot. The thin trace 
is the power spectrum, the thick dark line is the adopted model for this spectrum. The low-order polynomial modelling accounts for 
the confusion between red noise in the timing residuals and signal leakage caused by irregular sampling. The power spectra of PSRs 
J0437— 4715 and J1713+0747 have been modelled with quartics, while the spectrum of PSR J1857-I-0943 has been modelled with a cubic. 
The frequency coverage of PSR J0437— 4715 extends to much higher frequencies than those shown because of the very large number of 
timing residuals for this pulsar. 

data, we do not add sinusoids to the measured timing resid- 
uals and so a difTerent model for the power spectrum may 
be used because the spectral features are difTerent. 



